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Abstract 

A concurrent multiscale modeling methodology that embeds a molecular 
dynamics (MD) region within a finite element (FEM) domain has been 
enhanced. The concurrent MD-FEM coupling methodology uses statistical 
averaging of the deformation of the atomistic MD domain to provide 
interface displacement boundary conditions to the surrounding continuum 
FEM region, which, in turn, generates interface reaction forces that are 
applied as piecewise constant traction boundary conditions to the MD 
domain. The enhancement is based on the addition of molecular dynamics- 
based cohesive zone model (CZM) elements near the MD-FEM interface. 
The CZM elements are a continuum interpretation of the traction- 
displacement relationships taken from MD simulations using Cohesive Zone 
Volume Elements (CZVE). The addition of CZM elements to the concurrent 
MD-FEM analysis provides a consistent set of atomistically-based cohesive 
properties within the finite element region near the growing crack. Another 
set of CZVEs are then used to extract revised CZM relationships from the 
enhanced ESCM simulation of an edge crack under uniaxial loading. 
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I. Introduction 


Physics-based modeling of fracture begins at nanometer length scales in which atomistic 
simulation is used to predict the fonnation, propagation, and interaction of fundamental damage 
mechanisms. These mechanisms include dislocation fonnation and interaction, interstitial void 
formation, and atomic diffusion. The development of these damage mechanisms progress into 
microscale processes such as local plasticity and small crack fonnation. Ultimately, damage 
progression leads to macroscopic failure modes such as plastic yielding and large cracks 
exhibiting mixed-mode deformation. 

Modeling atomistic processes quickly becomes computationally intractable as the system size 
increases. With cunent computer technology, the computational demands of modeling suitable 
domain sizes (on the order of hundreds of atoms for quantum mechanics-based methods, and 
potentially billions of atoms for classical mechanics-based methods) and integrating the 
governing equations of state over sufficiently long time intervals quickly reaches an upper bound 
for practical analyses. In contrast, continuum mechanics methods such as the finite element 
method (FEM) provide an economical numerical representation of material behavior at length 
scales in which continuum assumptions are appropriate. Thus, multiscale modeling strategies are 
needed using the most efficient material model at different length scales for the purpose of 
developing a unified description of the hierarchy of processes that govern fracture. 

The concept of bridging length scales by coupling different computational paradigms is 
particularly attractive as a highly efficient means of reducing the computational cost of the 
simulations in cases that require the modeling of relatively large material domains to capture the 
complete deformation field, but where atomic and subatomic refinement is needed only in very 
localized regions (e.g., near a crack tip or dislocation core). Such computational issues arise in 
modeling crack nucleation and propagation, and in modeling dislocation formation and 
interaction. Coupled models allow the size limitations of the atomistic simulation to be avoided 
by embedding an inner domain, where complex dynamic processes and large deformation 
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gradients exist, within an outer domain, where the deformation gradients are small so that a 
continuum representation of the material becomes appropriate. 

Over the past decade, various methods that directly couple material representations at 
different levels of refinement have been developed and offer significant computational 
advantages compared to full atomistic simulations for predicting deformation and fracture 
processes. 1 ' 11 A common feature of all these approaches for coupling atomistic regions with 
continuum domains is the refinement of the FEM mesh down to the atomic scale However, a 
robust connection between continuum and discrete quantities can be achieved through a 
statistical averaging over scales where the discreteness of the atomic structure can be 
successfully approximated as a continuum. With the use of statistical averages to couple the two 
computational schemes, the developed approach constitutes a statistical coupling approach, and 
the developed MD-FEM coupling method is referred to as the embedded statistical coupling 
method (ESCM). ’ The resulting model that consists of an MD system embedded in a FEM 
domain is depicted in Figure 1. 

The statistical averaging over both time and volume of atomistic subdomains at the MD-FEM 
interface provides nodal displacement boundary conditions to the continuum FEM domain, 
which, in turn, generates interface reaction forces that are applied as constant traction boundary 
conditions 14 to the atoms within the localized MD subdomain. This process is iterated until the 
averaged quantities reach a stationary value. Typically, one finite element at the interface 
encompasses a region of several hundred to several thousand atoms. At this scale, the 
discreteness of the atomic structure is homogenized enough so that the FEM domain responds to 
the atomistic domain as an extension of the continuum. However, a moving crack that grows 
near the transition between the atomistic domain and the finite element region experiences an 
artificial constraint because using standard continuum finite elements along the crack path can 
not simulate the fracture process within the continuum. 

One means of accommodating the deformation within the finite element domain near the MD- 
FE interface is through the incorporation of cohesive zone model (CZM) elements in the 
analysis. Cohesive zone models approximate debonding processes along an interface 15 and are 
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frequently used in conjunction with finite element approaches to study fracture in a wide variety 
of materials. Recently, Cohesive Zone Volume Elements (CZVEs) 16 have been introduced in the 
MD analysis to calculate atomistic-based parameters defining CZMs. The resulting CZMs 
exhibit the aggregate behavior of a volume of atoms in the MD simulation. 16 ' 18 By placing the 
CZM elements along the plane of crack growth within the finite element region of the ESCM 
model, the fracture process can be simulated in the continuum as the crack approaches the MD- 
FE interface. Figure 2 depicts the coupled MD-FEM model containing an edge crack. CZM 
elements are identified that represent the cohesive properties of the GB in the continuum. Crystal 
I and II identify the MD regions that have been rotated through a tilt angle to form the GB. 

The present analysis considers intergranular fracture wherein MD simulations of fracture 
along grain boundaries (GBs) are used to obtain the atomistic response of crack propagation 
under applied loads along a characteristic length of the GB. Figure 3 shows a typical <1 10> 199 
symmetric tilt GB formed between two grains in a fee metal. Common neighbor analysis 
(CNA) 19 ' 20 is used to identify atoms in different crystallographic states: fee (small dots), hep 
(triangles), non-crystalline atoms (large dots), and surface atoms (squares), indicating existing 

vacancies in the GB. The length scale is in units of the lattice constant of Al, a 0 = 0.405 nm 

2 1 

using the potential of Mishin, et al. 

In this paper, an outline of the ESCM approach to coupling MD and FEM computational 
domains is presented in Section II. Section III summarizes the determination of CZM 
relationships from MD analysis using CZVEs to simulate fracture at higher length scales. Section 
IV presents the simulation of an edge crack propagating through a continuum into an atomistic 
domain where deformation along the crack plane in the continuum region is modeled using 
CZMs. This simulation provides atomistic resolution of near-tip processes in a crack with an 
extent that could not be simulated using purely MD methods. Finally, Section V presents a 
summary of the analysis and overall methodology. 
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II. The ESCM Model 


The structure of the basic ESCM model consists of four regions that are depicted in Figure 4. 
In this model, an Inner MD Region is connected to a FEM Domain along an interfacing segment 
in which the two material representations overlap. This segment is divided into an Interface MD 
Region that, together with the superposed FE nodes, constitutes the MD/FE Interface, and a 
Surface MD Region that does not interact with superposed FEM nodes but is used to compensate 
atomic free edge effects. The following subsections outline the function of each of these regions 
and the algorithmic coupling used in the ESCM. Complete details of the basic method are 

12 13 

presented by Saether et al. ’ 


II.l. MD and FEM model components 

The Inner MD Region is used to model material phenomena at the atomistic level. This Inner 
MD Region should be large enough to ensure a statistically smooth transition from an atomistic 
to continuum representation at the Interface MD Region while modeling any of the types of 
processes (e.g., mostly plastic deformation events such as dislocation or void nucleation, crack 
propagation) that evolve during the simulation. Together, the Inner, Interface and Surface MD 
Regions constitute the complete MD system. 

The only difference between the three MD regions is that, while the atoms in the Inner MD 
Region are subject only to the interatomic forces of neighboring atoms, the Interface and Surface 
MD Regions are also subjected to external forces involved in the ESCM coupling procedure. 

The fourth component is a FEM Domain that permits a large reduction in the computational 
cost of the performed simulations by replacing the atomistic representation with a continuum 
model in those regions of the system where the deformation gradients are small and atomic-level 
resolution is not necessary. 
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II. 2. The MD/FE Interface 


The first role of the MD/FE Interface is to provide a computational linkage between the MD 
region and FEM domain. The atoms that surround a given FEM node at the interface are 
partitioned to form a cell in the Interface MD Region, called an Interface Volume Cell (IVC), as 
shown in Figure 4. The IVCs are used to compute averaged MD displacements at the IVC mass 
centers that are then prescribed as displacement boundary conditions to the associated interface 
FEM nodes. During the coupled MD-FEM simulation, a spatial average within each k th IVC is 

performed to yield the center of mass displacement 8 ^ k , which is further averaged over a 
certain period of MMD time steps to yield the statistical displacement vector, 8f 4 ® 

■, M 

5" = («£*), =T7 £(W<;)-W «)) (1) 

M j= 1 

/ \ | N k / \ 

In the above expression, r CMk (/ / )= Xfi U / ) is the mass center of the k" IVC containing 

N k i'=t 

N/, atoms at positions r i at time tj of the / h MD step. The mass center displacement, 8qm k , in (1) 
is calculated relative to the initial zero-displacement position of the k 1b IVC, r CM k (0). 

The second role of the IVCs is to distribute reaction forces from the interface FEM nodes as 
external forces applied to the corresponding atoms within the IVC. The FEM reaction forces are 
computed as 


]M„)}+ [^]{M'„)}+P„]M,)} (2) 

where Kjy, Kjf, Kn are the terms in the partitioned FEM stiffness matrix representing the interior, 
far-field and interface regions, respectively. 8 V , 8/,- and 8 r are the corresponding displacement 
vectors with 8/ obtained from the MD region (Equation 1). 
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A third role of the IVCs, in the Interface MD region, together with the SVCs in the Surface 
MD region, is to incorporate a dissipative damping mechanism for phonons propagating into the 
interface. 13 

II. 3. The Surface MD Region 

In order for the MD domain to deform freely in response to applied reaction forces, it is 
modeled using free surface boundary conditions. However, the existence of a free surface 
introduces several undesirable effects in the MD system. First, it creates surface tension forces 
that must be removed to avoid distorting the MD response. Second, because atoms that should 
he within the cutoff radius of the free surface atom’s potential are missing, the coordination 
number of the surface atoms is reduced and they are, thus, less strongly bonded to the 
surrounding atomic field than those within the interior. Under sufficiently large reaction forces, 
these atoms may be separated from the surface layer causing a spurious surface degradation 
within the MD domain. To mitigate these free surface effects and to stabilize the atoms in the 
Interface MD Region, a Surface MD Region is introduced as shown in Figure 4. 

While the Surface MD Region eliminates the free surface effects within the Inner MD 
Region, it also introduces an undesirable fictitious stiffness which elastically constrains the 
deformation of the Inner MD Region. The separate effects of surface tension and the fictitious 
stiffness cannot be computed independently. However, their combined effect may be defined as a 
resultant force, f s , which acts at the boundary between the Surface and Interface MD Region, 
and is given by the sum of two components expressed as 

/,=! + *■ ( 3 ) 

In Equation 2, | is the elastic reaction of the Surface MD Region under deformation, and f is 
the force that results from the surface tension as shown in Figure 4. 


7 

American Institute of Aeronautics and Astronautics 



A way to mitigate both the surface tension and the elastic response of the Surface MD 
Region is to estimate and compensate for the force,/,. In the ideal case, when/ is fully 

compensated, the Surface MD Region acts as if it possesses zero stiffness and experiences no 
surface tension, thereby mitigating spurious influences on the Inner MD Region. Subdividing the 

Surface MD region into a number of SVCs helps to follow the variations of / along the 

perimeter of the Interface MD Region. For convenience, the partitioning of SVCs can be made to 
follow the IVC partitioning of the Interface MD Region. The resultant force is then calculated 

individually for each SVC. To compensate for these effects, a counterforce,/, is computed 

along the IVC/SVC interface and then distributed over the atoms of each SVC in a similar 
manner as the nodal reaction forces are applied to the IVCs of the Interface MD Region. The 

calculation of the counterforce, / , is presented in Saether, et al. ’ 

11. 4. MD-FEM coupling 

The MD-FEM coupling in the ESCM is achieved through an iterative equilibration scheme 
between the MD system and the FEM Domain. In this scheme, iterations begin with 
displacements at the MD/FE Interface that are calculated as statistical averages over the atomic 
positions within each IVC and averaged over the time of the MD simulation. These average 
displacements are then imposed as displacement boundary conditions, {</}, on the FEM 
domain. The resulting FEM boundary value problem is then solved to recover new interface 
reaction forces, { Rj } , resulting from the applied interface displacements and any imposed far- 
II eld loading. The new interface reaction forces, { R , }, are then distributed to the atoms in the 
IVCs, thus defining new constant traction boundary conditions on the MD system. This MD- 
FEM iteration cycle repeats until a stable equilibrium of both displacements and forces between 
the atomistic and continuum material fields is established at the interface. Details on calculating 
the interface reaction forces have been presented elsewhere. 12 ' 13 

11. 5. Constraint at MD-FEM interface 
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The ESCM was shown in References 12 and 13 to have many desirable features, including the 
applicability to three-dimensional domains, relatively coarse discretization of the FEM domain, 
and the ability to consider arbitrary temperature states. Additionally, the methodology was 
verified using a series of increasingly complex models that included the simulation of a 
propagating crack within the MD domain. However, a propagating crack that grows near the 
transition between the atomistic domain and the finite element region experiences an artificial 
constraint because the finite elements along the crack path cannot properly simulate the fracture 
process within the continuum, nor can they account for the deformation field (i.e., dislocation 
emission) produced by the approaching crack. This constraint places an inherent limitation on 
the original formulation of the ESCM and restricts its applicability to modeling an extent of 
crack propagation for which damage processes are entirely contained within the MD region. 
This limitation can be resolved by adding MD-based cohesive zone models (CZM) within the 
FEM domain along the anticipated plane of crack growth as shown in Figure 2. Section III 
summarizes the development of the CZM as formulated in Reference 16. 


III. Development of Traction-Displacement Relationships 

Through the definition of CZVEs shown in Figure 5, the results from MD simulations have 
been recast to obtain averaged continuum traction-displacement relationships that represent the 
cohesive interactions along a characteristic length of material. Using a purely MD model of a 
center crack along the 199 grain boundary presented in Reference 16, the stress and opening 
profiles are extracted along the crack growing in the system as shown in Figure 6. The figure 
shows distinct traction, cf yy {x), and displacement, A, profiles for the two crack tips of the center 
crack. The asymmetric behavior results from the orientation of the crystal lattice of the grains 
forming the grain boundary. A curve fit of f ' in Figure 6 is shown as a close approximation to 
the dyy{x) values extracted from the MD for the brittle tip and are consistent with the 
predominantly brittle fracture seen in this direction (left to right in Figure 5). In contrast, plastic 
deformation processes, including twinning and dislocation emission, were shown to dominate 
deformation at the ductile crack tip growing in the opposite direction (right to left in Figure 5). 16 
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The rfyyix) and A(x) profiles in Figure 6 represent the traction and opening displacement at a 
single instant of time during crack propagation. If the entire simulation of the interface 
debonding is divided into N t equal intervals of time t q C/e [0, N t ]), numerous similar profiles can 
be taken of many CZVEs placed along the GB. When plotted as a cfyy vs. A curve, each 
(cfyyixpjg), A (x p ,tq)) pair represents a point cj vy (k(x lh t q )). After sorting these data points in order 
of ascending A, so that A, < A i+ i, and taking a moving average (or a consecutive mean), 




i 

(2 M + 1) 


m , . 

2>MG+y) 

j=-M 


( 4 ) 


in which the results are averaged over M points backward and M points forward from A„ a 
construction of a statistically representative cohesive zone relationship z(A) can be made. 
Complete details of the modeling, fracture processes, and extraction of the CZM response are 
given in Yamakov, et al. 16 


IV. Crack Propagation through a Continuum-Atomistic Interface 

A continuum-atomistic model for crack propagation along the 199 GB shown in Figure 2 
and analyzed using the enhanced ESCM analysis is presented. The crack grows as an edge crack 
in the continuum environment subject to mode I loading conditions and penetrates into the 
atomistic environment (Figure 7). The continuity of the crack propagation process through the 
MD-FEM interface is achieved by using the CZM constitutive relation derived for interface 
debonding based on the pure molecular-dynamics simulations 16 and discussed in Section III. 
Note that the edge crack in Figure 7 is shown rotated to indicate that the simulation can be used 
for interfaces of arbitrary orientation. 

The dimensions of the coupled MD-FEM model used for this simulation are d fe = 900 nm in 
the x- and y- directions, with a circular MD domain in the x-y plane of diameter, duo = 45 nm. 
The system thickness in the z-direction is h = 2.9 nm with periodic boundary conditions to 

emulate a bulk atomic state in this direction. The simulation was performed at a temperature T = 

22 

300K and maintained by a Nose-Hoover thermostat. 
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An edge crack was inserted along the line of symmetry, c, in the FEM mesh which ended 
with four finite elements on the left side before reaching the MD-FEM interface (Figure 2). 
These four remaining FEs were connected through CZM elements placed along the crack 
direction (line c in Figure 2), which meet the atomistic GB at the MD-FEM interface. On the 
right side of the MD domain, the debonding interface is continued through another line of CZM 
elements, which extends to the opposite outer boundary of the FEM domain. 

The system was subject to uniaxial strain, £ yy , applied as displacement boundary conditions at 
the outer boundary of the FEM domain. The strain was gradually increased from 0.4% to 2% in 
steps of 0.4% every 20 MD-FEM iterations (20 ps MD simulation time), allowing the MD 
system to more smoothly follow the stress increase in the surrounding FEM domain. Figure 8 
presents the system configuration and <j vv stress state at the end of the second iteration when % = 
0.8%. 


The process of crack propagation is depicted in Figure 9 through a set of snapshots, 
magnified to show the physical configuration of the MD domain in the center of the MD-FEM 
coupled model. The GB interface in the MD domain is visualized as a line of blue and red atoms, 
indicating atoms in a disordered or hep surrounding, respectively, reflecting the non-crystalline 
structure of the GB interface. In Figures 9(a) and 9(b), the crack propagates through the FEM 
mesh, sequentially opening the CZM elements placed on its path. In Figure 9(c), the crack tip 
enters the MD domain. The snapshot has captured the instant when the debonding of the Surface 
and the Interface MD Regions is accomplished in unison with the opening of the last CZM 
element. The crack tip is about to enter the inner MD region and to continue its propagation in a 
fully atomistic way. In Figure 9(d), decohesion of the GB in the inner MD region has begun. 
Shortly afterward, in Figure 9(e), plastic processes start to develop in terms of nucleation of a 
pair of dislocations, seen as two short blue lines, marking the extended dislocation cores on both 
sides of the crack tip. Figure 9(f) presents the crack propagated well inside the MD domain. 
Active dislocation nucleation has already developed and accompanies the crack’s propagation. 
An important result of this process, shown in Figures 9(a-f), is the preservation of the structural 
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continuity between the atomistic and the continuum domains at the MD-FE interface through the 
entire process of crack propagation from the FEM to the MD systems. 

The stress continuity along the MD-FE interface is also well preserved during the entire 
process. Figure 10(a-f) shows the stress state corresponding to the configurational snapshots 
presented in Figure 9. The largest discrepancy between the stress of the MD system and the FEM 
mesh is seen in Figures 10(b) and 10(c) where the dynamics of the MD domain did not converge 
completely before the next stepwise increase of the strain in the static FEM domain was applied. 
The continuity is improved later in the simulation, in Figures 10(d-e), when the crack propagates 
further within the MD region. 

The problem of coupling the dynamic response of the MD domain to the static response of 

n 

the FEM system is discussed in detail in Saether, et al. ~ While the dynamic coupling in the 
ESCM approach cannot be fully achieved without the use of a dynamic FEM simulation, the 
iterative scheme of coupling, in which the FEM state is continuously updated in accordance with 
the reported MD displacements, appears to be sufficient to reproduce the evolution of systems 
with relatively slow dynamics. 

In the original purely MD simulation 16 that was used to obtain CZM relationships, model size 
restrictions resulted in a finite geometry distortion. This effect was alleviated using the enhanced 
ESCM simulation, and the brittle Mode I CZM was recomputed. The procedure for the 
determination of the CZM parameters in a purely MD simulation is the same as that used in the 
coupled MD-FEM simulations. The insert in the bottom half of Figure 1 1 shows the schematic 
diagram of the discretization of the system volume in the MD simulation and the defining of the 
CZVEs for extracting the parameters for the CZM elements suitable for use in a continuum 
simulation. A contrast of the CZVE configurations for analyzing the embedded crack using a 
purely MD analysis and for the edge crack analyzed using ESCM is seen by comparing Figures 5 
and 1 1 . A contrast between the stress and opening profiles for both simulations can be seen by 
comparing Figures 6 and 11. 
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The use of the coupled MD-FEM ESCM approach eliminated the finite geometry distortion 
seen in the purely MD simulation. 16 A comparison of the recomputed brittle Mode I CZM with 
the initial pure MD results from Reference 16 is presented in Figure 12. As shown, the peak 
stress is increased slightly but the overall energy and functional form of the CZM is closely 
reproduced. 


V. SUMMARY 

A coupled MD-FEM analysis of intergranular fracture based on the embedded statistical 
coupling method (ESCM) is presented. The ESCM approach is based on solving the boundary 
value problem through an iterative procedure for both MD and FEM systems at their common 
interface. The FEM system is loaded along the MD-FEM interface by nodal displacement 
boundary conditions obtained as statistical averages of the atomic positions in the MD system at 
the mass centers of associated Interface Volume Cells (IVCs). The MD system is simulated 
under constant traction boundary conditions obtained from the FEM system as reaction forces to 
the MD displacements at the interface. This approach allows the continuity at the MD-FEM 
interface to be achieved at different length scales inherent to both systems. 

An application of ESCM is presented that utilizes Cohesive Zone Volume Elements (CZVEs) 
and a moving average approach to determine CZM traction-displacement parameters of a 199 
GB in aluminum. The model incorporated an edge crack with the tip propagating into the Inner 
MD Region of ESCM. The continuity of the crack propagation process through the MD-FEM 
interface was achieved by the derivation of the CZM constitutive relation for interface debonding 
using CZVEs and moving averages based on an earlier pure molecular-dynamics simulation. 
CZM parameters for Mode I fracture were determined and compared with those determined from 
the original pure MD simulation. For a ‘brittle’ crack tip, the peak stress was increased slightly in 
the ESCM simulation but the overall energy and functional form of the CZM was closely 
reproduced. 

The presented technique demonstrates the versatility of the ESCM approach and opens new 
opportunities for performing multiscale atomistic-continuum modeling of fracture processes 
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where the fracture region is not limited to the atomistic domain. Additionally, casting atomistic- 
level failure processes into continuum constitutive relations is a necessary step towards 
developing a truly physics-based multiscale methodology to extend material simulation from 
nanoscopic to microscopic length scales. 
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Figure 1. An embedded MD region within a FEM domain. 
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Figure 2. Configuration of the simulation model for edge-crack propagation 
along a grain-boundary in aluminum. 
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Figure 3. Crystallography and structure of the GB interface. 
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Figure 4. Structure of the coupled MD-FEM ESCM approach showing the forces 
along the interface. 
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Figure 5. Schematic diagram of the system volume slicing in the MD 
simulation of a 199 grain boundary and defining the cohesive zone volume 
elements (CZVE) for extracting cohesive-zone model (CZM) parameters. 
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Figure 6. Stress and opening profiles extracted along the embedded crack growing in a 
purely MD simulation. A snapshot of the corresponding crack configuration is shown 
at the bottom. 
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Figure 7. CZVE elements along the crack path inside the inner MD region 
within the ESCM. 
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Figure 8. ESCM simulation of an edge crack propagating along a <1 10>/E99 
symmetric tilt grain boundary in aluminum for 2% far-field uniaxial strain. 
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Figure 9. Snapshots of the ESCM model configurations simulating an edge crack propagating 
along a GB in a bicrystal of aluminum with CZM elements along the GB interface. 
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Figure 10. Stress map of <jyy corresponding to the snapshots given in Figure 9. 
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Figure 11. Stress and opening profiles extracted from cohesive-zone-volume-elements 
(CZVE) placed along the path of an edge crack growing along a grain-boundary inside 
the molecular dynamics region. 
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Figure 12. Comparison between the brittle CZM derived using 
pure MD and ESCM models. 
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